function [HHstates,wage_moments,transitions]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,wi_1,Wf,Wi,Wn)


df_1=Par1(1);
di_1=Par1(2);
lnf_1=Par1(3);
lni_1=Par1(4);
lfi_1=Par1(5);
lif_1=Par1(6);
alphaf=Par1(15);
betaf=Par1(16);
alphai=Par1(17);
betai=Par1(18);
% Data simulation (MSM)

N_households=10000; % # households
maxtime=400;   % # quarters
rng('default'); % set seed

% Pre-allocation
value=zeros(N_households,maxtime);
state=zeros(N_households,maxtime);
wage_draw=ones(N_households,maxtime);

for iter=1:N_households
    
    time=1; % households are born in the situation 1: "N"
    scale = lni_1  + lnf_1  + (1 - lni_1)*(1-lnf_1);
    randnumber=scale*rand;
    
    if randnumber < lnf_1
        urn = random('Beta', alphaf, betaf);
        transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
        wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
        value(iter,time)=Wf(wage_aux)*(Wf(wage_aux)>Wn)+Wn*(Wf(wage_aux)<=Wn);
        state(iter,time)=(Wf(wage_aux)<=Wn)+2*(Wf(wage_aux)>Wn);
        wage_draw(iter,time) = wage_aux*(Wf(wage_aux)>Wn) + (Wf(wage_aux)<=Wn);
        
    elseif randnumber >=lnf_1 && randnumber<(lni_1+lnf_1)
        urn = random ('Beta', alphai, betai);
        transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
        wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
        value(iter,time)=Wi(wage_aux)*(Wi(wage_aux)>Wn)+Wn*(Wi(wage_aux)<=Wn);
        state(iter,time)=3*(Wi(wage_aux)>Wn)+ (Wi(wage_aux)<=Wn);
        wage_draw(iter,time) = wage_aux*(Wi(wage_aux)>Wn) + (Wi(wage_aux)<=Wn);
         
    else
        wage_draw(iter,time)=1;
        value(iter,time)=Wn;
        state(iter,time)=1;
        
    end
    
    
    time=2;
    
    while time<=maxtime
        
        
        
        
        %% F in t-1
        if (state(iter,time-1) == 2)
            scale = df_1 + lfi_1 + (1 - df_1)*(1 - lfi_1);
            randnumber=scale*rand;
            if  (randnumber < df_1)
                wage_draw(iter,time)=1;
                value(iter,time)=Wn;
                state(iter,time)=1;
            end
            if randnumber >= df_1 && randnumber<(lfi_1+df_1)

                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                value(iter,time)=Wi(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=3*(value(iter,time)> value(iter,time-1))+2*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw(iter,time) = wage_draw(iter,time - 1);
                else
                    wage_draw(iter,time) = wage_aux;
                end
            end
            if randnumber>= (lfi_1+df_1)
  		wage_draw(iter,time)=wage_draw(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
            end
        end
        %% I em t-1
        if state(iter,time-1)==3
            scale = di_1 + lif_1 + (1-di_1)*(1-lif_1);
            randnumber=scale*rand;
            
            if randnumber < di_1
                wage_draw(iter,time)=1;
                value(iter,time)=Wn;
                state(iter,time)=1;
            end
   
            if randnumber >= (di_1) && randnumber<(lif_1+di_1)
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                value(iter,time)=Wf(wage_aux);
                value(iter,time)=value(iter,time)*(value(iter,time)>value(iter,time-1))+value(iter,time-1)*(value(iter,time)<=value(iter,time-1));
                state(iter,time)=2*(value(iter,time)> value(iter,time-1))+3*(value(iter,time)<= value(iter,time-1));
                if state (iter,time) == state(iter,time-1)
                    wage_draw(iter,time) = wage_draw(iter,time - 1);
                else
                    wage_draw(iter,time) = wage_aux;
                end
            end
            if randnumber>=(lif_1+di_1)
                wage_draw(iter,time)=wage_draw(iter,time-1);
                value(iter,time)=value(iter,time-1);
                state(iter,time)=state(iter,time-1);
            end
        end
       

        % N in t-1
        if state(iter,time-1) ==1
            scale = lnf_1 + lni_1 +(1 - lni_1)*(1 - lnf_1);
            randnumber=scale*rand;
            
            if randnumber < lnf_1
                urn = random('Beta', alphaf, betaf);
                transf_wage = wf_1(1) + urn*(wf_1(50)-wf_1(1));
                wage_aux=find(abs(wf_1 - transf_wage)<=min(abs(wf_1 - transf_wage)));
                value(iter,time)=Wf(wage_aux)*(Wf(wage_aux)>Wn)+Wn*(Wf(wage_aux)<=Wn);
                state(iter,time)=2*(Wf(wage_aux)>Wn)+(Wf(wage_aux)<=Wn);
                 if state (iter,time) == state(iter,time-1)
                    wage_draw(iter,time) = wage_draw(iter,time - 1);
                else
                    wage_draw(iter,time) = wage_aux;
                end
            end
            if randnumber >=lnf_1 && randnumber<(lni_1+lnf_1)
                urn = random ('Beta', alphai, betai);
                transf_wage = wi_1(1) + urn*(wi_1(50)-wi_1(1));
                wage_aux=find(abs(wi_1 - transf_wage)<=min(abs(wi_1 - transf_wage)));
                value(iter,time)=Wi(wage_aux)*(Wi(wage_aux)>Wn)+Wn*(Wi(wage_aux)<=Wn);
                state(iter,time)=3*(Wi(wage_aux)>Wn)+ (Wi(wage_aux)<=Wn);
                 if state (iter,time) == state(iter,time-1)
                    wage_draw(iter,time) = wage_draw(iter,time - 1);
                else
                    wage_draw(iter,time) = wage_aux;
                end
            end
            
            if randnumber>=(lni_1+lnf_1)
                value(iter,time)=Wn;
                state(iter,time)=1;
            end
        end
        
        
        time=time+1;
        
    end
    
end

% Moments calculation

% use only last 2 quarters of simulation for each household
NF=sum(state(:,end-1)==2)/N_households;
NI=sum(state(:,end-1)==3)/N_households;
NN=sum(state(:,end-1)==1)/N_households;


HHstates=[NF;NI;NN]; % joint-job state

wage_formal_1=wf_1(wage_draw(:,end-1));
wage_informal_1=wi_1(wage_draw(:,end-1));


indicewf_1=find(state(:,end-1)==2);
wage_formal_1=sort(wage_formal_1(indicewf_1));
indicewi_1=find(state(:,end-1)==3);
wage_informal_1=sort(wage_informal_1(indicewi_1));

sdwi_1 = std(wage_informal_1);

sdwf_1 = std(wage_formal_1);


meanwf_1=mean(wage_formal_1);

p10wf_1=prctile(wage_formal_1,10);
p25wf_1=prctile(wage_formal_1,25);
p50wf_1=prctile(wage_formal_1,50);
p75wf_1=prctile(wage_formal_1,75);
p90wf_1=prctile(wage_formal_1,90);


meanwi_1=mean(wage_informal_1);

p10wi_1=prctile(wage_informal_1,10);
p25wi_1=prctile(wage_informal_1,25);
p50wi_1=prctile(wage_informal_1,50);
p75wi_1=prctile(wage_informal_1,75);
p90wi_1=prctile(wage_informal_1,90);



wage_moments=[meanwi_1; meanwf_1;  sdwi_1; sdwf_1;  p10wi_1 ; p10wf_1; ...
    p25wi_1; p25wf_1;  p50wi_1; p50wf_1;  p75wi_1; p75wf_1;...
    p90wi_1; p90wf_1]; % wages

% conditional transitions
formal_1_time1=(state(:,end-1)==2);
informal_1_time1=(state(:,end-1)==3);
nonemp_1_time1=(state(:,end-1)==1);


formal_1_time2=(state(:,end)==2);
informal_1_time2=(state(:,end)==3);
nonemp_1_time2=(state(:,end)==1);


Dfn_1=sum(formal_1_time1==1 & nonemp_1_time2==1)/sum(formal_1_time1==1);
Din_1=sum(informal_1_time1==1 & nonemp_1_time2==1)/sum(informal_1_time1==1);
Dnf_1=sum(nonemp_1_time1==1 & formal_1_time2==1)/sum(nonemp_1_time1==1);
Dni_1=sum(nonemp_1_time1==1 & informal_1_time2==1)/sum(nonemp_1_time1==1);
Dfi_1=sum(formal_1_time1==1 & informal_1_time2==1)/sum(formal_1_time1==1);
Dif_1=sum(informal_1_time1==1 & formal_1_time2==1)/sum(informal_1_time1==1);

transitions=[Dfn_1;Din_1;Dnf_1;Dni_1;Dfi_1;Dif_1];
end





